Persistent current formation in a high-temperature Bose-Einstein condensate: 
an experimental test for c-field theory 
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We quantitatively model a stirred Bose-Einstein condensate at high temperature using grand-canonical c- 
field theory, comparing simulations with available experimental data. Stirring generates a chaotic array of 
quantum vortices that decays, via thermal damping and noise, to form a macroscopic persistent current in a 
toroidal confining potential. We perform 3D numerical simulations of the full dynamical evolution of the system 
within the truncated- Wigner representation using a consistent energy cutoff and a priori determined reservoir 
parameters. This work gives the first quantitative comparison of theory with experiment for dissipative vortex 
dynamics in a trapped Bose-Einstein condensate from first principles, with no fitted parameters. We find that 
both damping and noise are required to give a quantitative account of the dynamics, and that the size of the 
persistent current and timescale of its formation are in close agreement with experimental data. In contrast, 
comparison of experimental data with damped Gross-Pitaevskii theory excludes it as a quantitative model of 
high temperature vortex dynamics. 
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As highly controllable degenerate matter wave systems, di- 
lute Bose gases offer a unique window into the realm of many 
body quantum mechanics Q. Yet developing a quantitative 
non-equilibrium description of high-temperature Bose gases 
poses a major theoretical challenge GIO, particularly near the 
critical temperature where the breakdown of mean field theory 
renders two-fluid theories (4) inoperative. Exact methods of- 
fer insight for small systems [5 ], but for large hot systems, that 
typically occur near the critical point of evaporative cooling, 
quantitative tests of theory have been confined to equilibrium 
critical fluctuations [6], and collective modes (7). As conspic- 
uous experimental signatures of superfluid motion with long 
lifetimes |8], quantum vortices offer a different test for non- 
equilibrium theory. 

Since the first observations of quantized vortices in dilute 
gas Bose-Einstein condensates (BEC) [9 ], experimental stud- 
ies of their properties have proliferated ifTOl . Vortex motion is 
highly sensitive to thermal fluctuations and nonlinear interac- 
tions between vortices and other excitations in the fluid ifTTIl . 
Many experiments depend on the presence fl2l or forma- 
tion fT3l of a rotating thermal reservoir to take the superfluid 
from a non-rotating state to one containing vortices. Dissi- 
pative vortex dynamics in harmonically trapped BECs have 
been treated within Zaremba-Nikuni-Griffin (ZNG) theory at 
low temperatures |[T4l . and the Stochastic Projected Gross- 
Pitaevskii equation (SPGPE) [ 15 ] at high temperatures l8lfT6l, 
but the predictions for vortex lifetimes have yet to be tested 
experimentally. In an annular trapping geometry, such as may 
be created using an obstacle potential in a harmonic trap ifTTlL 
a vortex can be freely moving in the bulk superfluid, or pinned 
to the obstacle, forming a stable persistent current. Persistent 
currents have previously been created in BECs through coher- 
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ent optical manipulation 11181 , and the Kibble Zurek mecha- 
nism [17]. Despite the importance of quantized vortices in 
in dissipative superfluid phenomena and in the formation and 
stability of persistent currents, their motion at high tempera- 
ture is largely untested. 

In this article we present a study of persistent current for- 
mation via dissipative vortex dynamics in the presence of a 
large thermal reservoir. We perform large scale numerical 
simulations of an experimental forcing sequence |19] using 
the SPGPE to provide the first quantitative test of theory of 
vortex dynamics at high temperature, with no fitted param- 
eters. We also compare the SPGPE with the damped GPE 
(dGPE) l20l which is obtained by setting the noise in the 
SPGPE to zero, and has been a useful tool for understanding 
vortex motion lf2T1l22ll . 

The SPGPE is a grand-canonical c-field method formulated 
from a microscopic derivation of reservoir interactions in the 
Wigner phase-space representation (15] E3, and is valid right 
through the phase transition [3 ]. The equation of motion bears 
a superficial resemblance to the GPE, but also contains damp- 
ing and noise terms arising from the reservoir interaction, and 
evolves both the condensed and non-condensed fractions of 
the system lying below a specified energy cutoff. The cutoff 
is a central formal and technical aspect of the theory. It allows 
the formally divergent continuum SGPE theory to be consis- 
tently extended beyond one spatial dimension fi31 . In this 
work we use the simple-growth SPGPE (24), which neglects 
reservoir interactions involving number-conserving scattering 
between c-field and thermal reservoir atoms. These processes 
are known to be weak in quasi-equilibrium situations 0. In 
a simple-growth SPGPE study of vortices occurring sponta- 
neously during the phase-transition ifTvTl a single fitted param- 
eter (the reservoir coupling strength) was used to give vortex 
formation data in close agreement with the experiment. How- 
ever, in near-equilibrium situations the reservoir interaction 
parameters can be determined a-priori (SEU, allowing the 
SPGPE to perform quantitatively accurate calculations of dis- 
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FIG. 1. Schematic of the experimental sequence and simulation pa- 
rameters, (a) The initial state consists of a trapped Bose gas at 
T ~ 0.9Tc in a cylindrically- symmetric harmonic trap augmented 
with an optically induced Gaussian obstacle potential. The potential 
is constant in the z direction and coincident with the symmetry axis 
of the harmonic trap. In- situ absorption images of the experimental 
atom density are shown in the transverse (b) and axial (c) directions. 
(d) The center of the harmonic trap is induced to complete a single 
revolution around the z-axis. The obstacle beam executes one cir- 
cular orbit centered at (x(t),y(t)) — ro(l — cos «£, sin at), with 
ro = 2.875/xm, k = 6ty s _1 . This is modeled in the frame of the 
obstacle beam, with an incoherent (I-) region coupled to a coher- 
ent (C-) region described with the SPGPE (shown schematically, see 
text), (e) The dissipation is a significant effect for long time-scales 
(during the experimental hold). Absorption images after significant 
hold times show either (/) many vortex cores after an additional ra- 
dial expansion stage, or (g) a large density minimum corresponding 
to a persistent current (without radial expansion). 



sipative Bose-gas dynamics at high-temperature, provided the 
incoherent field is not significantly disturbed. 

In our c-field description the truncated- Winger field is ex- 
panded on a basis of harmonic oscillator eigenstates. The 
modes of the system are divided into two distinct regions: the 
coherent region (C) consisting of modes with energy less than 
a specified cutoff (e cut ), and the incoherent region (I) which 
contains the remaining high energy, quasi-equilibrium states. 
The I-region acts as a thermal reservoir for the C-region, and 
is assumed to be in thermal equilibrium at a temperature T 
and chemical potential /i and is described by a semiclassical 
Bose-Einstein distribution. The C-region is treated using the 
truncated Wigner method, where accounting for interactions 
between the C- and I-regions leads to a stochastic differential 
equation for the C-region dynamics. The equation of motion 
is 1241 



(1) 



#c(r,t) = V c { - l -L c ^ c (v,t)dt 
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(/x-L c )^c(r,t)dt + d^(r,t)}, 



where the projection operator Vc implements the en- 
ergy cutoff in the spectral basis, and the complex 
Gaussian noise satisfies (dW(r, t)dW(r', t)) = 0, 
(dW*(r,t)dW(r',t)) = (2jk B T/h)5 c (r,r')dt, where 
£c(r,r') = S n ec^( r )0n( r/ ) is a delta-function in C. 
Lq is the Hamiltonian evolution operator for the C-region 
Lc^c = (#sp + d\^c\ 2 ) ipc, where the single-particle 
Hamiltonian is H sp = -h 2 V 2 /2m + V(r, t), g = 47rh 2 a/m 
characterizes the strength of the atomic interaction, and a is 
the s-wave scattering length. 



The first term of the SPGPE ([T]) describes Gross-Pitaevskii 
evolution of the C-region, while the second and third terms 
account the growth of the C-region from scattering of two I- 
region atoms, and the corresponding time reversed process. 
Experiments are most commonly described by the system 
temperature and total atom number Nt. In the SPGPE im- 
plementation, the choice of /i(T, Nt) controls the total atom 
number Nt, while the choice of e cut (T, Nt) dictates the oc- 
cupation at the cutoff which must be of order unity for the 
classical field description of the C-region to be valid 0. 
A Hatree-Fock method can be used to accurately estimate 
SPGPE parameters for harmonically trapped systems close to 
equilibrium O . The growth rate 7 can be calculated for near 
equilibrium situations where the incoherent region is well de- 
scribed by an ideal semiclassical Bose-Einstein distribution, 
where it takes the form 1 24 1 



7(T,/i,e cut ) = 70 



00 e /3/i(fc+l) 



2/3e cu t k 
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p/3e C ut 



,1,* 



(2) 



where $[u,v, w] is the Lerch transcendent and the dimension- 
less rate constant is 70 = 8a 2 /\ 2 dB for de Broglie wavelength 
XdB — \/27rh 2 /rnksT. This explicit calculation of 7 means 
all SPGPE parameters are found in a physically consistent 
manner prior to simulation, giving a first-principles treatment 
of damping with no fitted parameters. 

We describe the toroidal system using a harmonic- 
Gaussian external potential V(r,t) = Vho{ y ) + Vb(r,£), 
where V H o(?) = f [ul(x 2 + y 2 ) + ou 2 z 2 ] is the 
harmonic oscillator potential, and the time depen- 
dent Gaussian potential is given by Vcf(r,t) = 
Vbexp [-[(x - x(t)) 2 + (y - y{t)) 2 )/ a 2 ]. To match 
the persistent current formation experiment we have 
(uJr^z) = 2tt x (8,90) Hz, a = 23 /xra/>/2, and 
Vo = 5ShCd, where Co = (cj^cj^) 1 / 3 . Note this poten- 
tial has a minimal effect on the I-region as it is well 
contained within the C-region, since cfq <C i? C ut> where 
R cut = yj 2e cut / muj 2 = 73 /im is the spatial cutoff im- 
posed by e cut . Thus our chosen basis gives a complete 
representation of the C-region. 

To model the persistent current experiment we create an 
initial equilibrium state by evolving Eq. ([T]) with parame- 
ters (T, /i, e cut ) = (98nK, 34fkD, 83hCu) with the Gaussian 
potential at (x,y) = (0,0). These self consistently deter- 
mined parameters sample the equilibrium ensemble of 
the SPGPE, giving Nt = 2.6 x 10 6 total atoms (23, to 
match the experimental value. For dynamics Eq. ^ gives 
7(98nK, 34to, 83hCu) = 8 x 10~ 4 and we model the stir by 
setting (x(t),y(t)) = r*o(l — cos(ftt), sin(/ct)), to move the 
Gaussian beam in a single circle of radius ro = 2.875 /im, 
about the point (x,y) = (7*0, 0), with angular frequency 
k = 27r/(333 ms) = 6tt s _1 (see Fig. [I]). During the 
post-stir hold period where dissipation is a significant effect 
[Fig.[T](e)], there is an additional cooling stage of the exper- 
iment |19l . Modeling of this process via a time dependent 
reservoir revealed no significant modification to the dissipa- 
tive evolution, and thus the details of cooling are omitted here. 
Our numerical method is a spectral Galerkin method based on 
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FIG. 2. Column densities and phase slices (through the z — plane) 
showing the c-field dynamics of modeling the persistent current for- 
mation experiment, (a) DPGPE evolution, where dissipation on its 
own leads to long vortex lifetimes, (b) SPGPE evolution for a sin- 
gle trajectory, where thermal fluctuations greatly increase the vortex 
decay rate and enhance the size of the persistent current. The red 
circles at t = s show the boundary of the region within which we 
count vortices. 



Gauss-Hermite quadrature [ 3 , 26 ]. Due to large particle num- 
ber, the trap oblateness, and the need for a high energy cutoff, 
we require ~ 10 5 modes in the C-region, making each trajec- 
tory numerically challenging [27]. Our implementation of the 
DPGPE is identical to that of the SPGPE, without the noise in 
Eq. Q. 

Fig. [2] shows a comparison of individual trajectories of the 
DPGPE [Fig. [2] (a)] and the SPGPE [Fig.[2](b)]. The dynam- 
ics of both methods are qualitatively similar. Multiple vortices 
are nucleated during the stirring procedure, which are most 
easily identified in the phase profile. These vortices evolve in 
a complicated manner, where they decay through a range of 
processes: decay to the exterior condensate boundary, inter- 
nal vortex/antivortex annihilation, or via pinning at the cen- 
tral potential. Finally the system evolves into a stable persis- 
tent current with no free vortices in the bulk fluid. However 
the quantitative differences between the DPGPE and SPGPE 
are readily apparent from the timescale of the decay of free 
vortices, and the size of the macroscopic persistent current 
formed. In Fig. [2] (a) we see the DPGPE requires more than 
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FIG. 3. Average number of vortices N v <[3j as a function of time t, 
where the inset shows the long time dynamics, th corresponds to 
the end of the 0.333 s stir procedure. The SPGPE ensemble average 
is shown for for V = 58to (SPGPE (a), red) , and V = 67to 
(SPGPE (b), blue), where in both cases a stable persistent current 
is formed after t ~ 3 s. The DPGPE simulation (gray) develops a 
stable persistent current after t ~ 20 s. The experimental data is 
shown by the squares with error bars. 



30 s to evolve into a stable persistent current with a circula- 
tion of 3 quanta. In contrast, in Fig.[2](b), SPGPE evolution 
generates a stable charge-4 persistent current after 2.8 s. The 
full evolution of one trajectory of the SPGPE is shown in Sup- 
plementary Movie SI 1281 . 

We compare the DPGPE and SPGPE by calculating the 
number of vortices at a given time in the numerical simula- 
tions (in SPGPE we average over 16 trajectories), and compar- 
ing this with the number found in the experiment. In counting 
vortices, we limit our region of detection 7Z to radii within 
22.9 fim < 1Z < 34.4 jira centered about the center of the 
trap, to avoid counting thermal fluctuations. This region is 
shown in Fig. [2] by the red circles at t = s. Experimen- 
tally the vortices were counted after ramping down the Gaus- 
sian beam, so this number includes any vortices pinned to the 
Gaussian beam. Thus we calculate the number of vortices 
from the numerical simulations via 



N v = N f + N T = N- 



f 



<£ v(r) • d\ 



(3) 



where Nf is the number free of vortices in the region of de- 
tection, and the Bj is the inner boundary defining the num- 
ber pinned (iVj). Eq. ^ gives the number of vortices in 
the toroidal condensate including those pinned to the Gaus- 
sian potential, thus it can be compared with experimental ob- 
servations. In Fig. [3] we plot N v for the DPGPE simula- 
tion (gray curve), and for the SPGPE ensemble average (red 
curve). We have also included an extra set of SPGPE trajecto- 
ries (blue curve), where we have used a higher barrier height 
of Vo = 67 HO to account for the experimental uncertainty in 
this parameter. In this set of simulations we use SPGPE pa- 
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rameters of (T, /i, e cut , 7) = (98nK, 35 to, 84to, 8 x 10~ 4 ) 
to give the same atom number and temperature as for the sim- 
ulations with the lower barrier height of Vo = 58 to. Fig. [3] 
shows an initial rise in N v due to nucleation of vortices from 
the stir. More vortices are nucleated in the SPGPE simula- 
tions. The peak number of vortices occurs at t = 0.21 s in 
both SPGPE ensemble averages, after which there is a dra- 
matic drop in N v . In comparison, the peak number of vor- 
tices occurs at t = 0.33 s in the DPGPE, and N v decreases 
more gradually. We compare these calculations with available 
experimental data (black squares). Both SPGPE calculations 
agree well with the first experimental data point at t = 6.33 s, 
lying within the experimental uncertainty. At this stage N v 
has reached a stable value for the SPGPE simulations, where 
all free vortices have left the condensate and the system is in a 
stable persistent current. In contrast, there are still many free 
vortices after t = 6.33 s in the DPGPE simulation. 

The quantitative difference at long evolution times is shown 
in the inset of Fig. [3] where we compare the data of simula- 
tions and experiment. Experimentally, N v decays due to a 
loss of pinned vortices caused by a slow drift in the mag- 
netic trap center |fT9l . This means the drop in N v is due to 
the lack of complete control of the trap, rather than decay of 
free vortices. This is also supported by experimental column 
densities where we observe no evidence of free vortices for 
t > 10 s (29). Thus the agreement of N v (at t = 6.33 s) be- 
tween the SPGPE calculations and experiment show that the 
SPGPE accurately predicts the size of the persistent current, 
and the decay-time of free vortices in the condensate. Free 
vortices exist in the DPGPE simulation for the first 20 s of 
evolution, over which time N v is never comparable with ex- 
periment. After t = 22 s, the DPGPE stabilizes into a charge 3 
persistent current. Note that this is smaller than that seen with 



the SPGPE, and is only similar to the experimental persistent 
current for times longer than 20 s. 

The vortex detection region is limited by the ability to see 
free vortices in competition with thermal fluctuations in the 
SPGPE simulations. Thus the absence of thermal fluctuations 
in the DPGPE simulation means that N v found with our spec- 
ified 1Z is a lower bound. In Fig. [2] (a), we see there are free 
vortices present at t = 24 s, despite Fig. [3] showing that the 
system has reached a stable persistent current at this time. 

In summary, we have modeled the persistent current exper- 
iment of Neely et al fT9l using the grand-canonical c-field 
method. The SPGPE theory of reservoir interactions predicts 
a persistent current formation time and winding number that 
closely matches the experiment. Such quantitative accuracy 
is not provided by the damped GPE, and the regime of sig- 
nificant vortex dissipation in the experiment (T « 0.9T C ) is 
well outside the region of validity of generalized mean-field 
theories. In general, both damping and noise are required to 
give a quantitative description of dissipation in open quan- 
tum systems, and our results demonstrate the central impor- 
tance of thermal noise in high temperature Bose gas dynam- 
ics. We have demonstrated the crucial effect thermal fluctu- 
ations have on the physical decay of quantized vortices, and 
performed an a priori quantitative test of the SPGPE. Our ap- 
proach provides a general quantitative framework for studying 
high-temperature dynamics in three-dimensional Bose gases. 
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